In Part 1, we showed how to preprocess scRNA-seq isoform datasets and create a multi Assay Seurat object. In this part, we’ll use this Seurat object to perform differential isoform expression analysis and explore pertinent biological questions. We’ll focus on searching for genes with a transcript isoform expression that differed between clusters, which we refer to as “isoform switches”.
We start the analysis by loading the necessary packages:
library(dplyr)
library(patchwork)
library(ggplot2)
library(reactable)
Plotting functions
see_isoforms_expression = function(sobj, gene, limits) {
isoforms = grep(paste0("^", gene, "\\."), rownames(sobj), value = TRUE)
plot_list = lapply(isoforms, FUN = function(one_isoform) {
p = Seurat::FeaturePlot(sobj, features = one_isoform,
reduction = name2D, order = TRUE) +
ggplot2::labs(title = one_isoform) +
ggplot2::scale_color_gradientn(limits = limits, colors = c("lightgray", "#FDBB84", "#EF6548", "#7F0000", "black")) +
ggplot2::theme(aspect.ratio = 1,
plot.subtitle = element_text(hjust = 0.5)) +
Seurat::NoAxes()
return(p)
})
return(plot_list)
}
see_genes_expression <- function(gene){
Seurat::FeaturePlot(seurat_obj, gene, reduction = name2D)+
ggplot2::scale_color_gradientn(limits = c(0,4), colors = c("lightgray", "#FDBB84", "#EF6548", "#7F0000", "black")) +
ggplot2::theme(aspect.ratio = 1,
plot.subtitle = element_text(hjust = 0.5)) +
Seurat::NoAxes()
}
We set the global settings of the analysis. We will store data there :
out_dir = "./output/"
data_dir <- "./data/"
We load the Seurat object:
seurat_obj = readRDS(paste0(out_dir, "/datasets/", "seurat_obj_annotated.rds"))
color_markers = readRDS(paste0(out_dir, "/datasets/", "color_markers.rds"))
This is the projection of interest:
name2D = "RNA_pca_10_tsne"
We visualize cells:
cell_type = Seurat::DimPlot(seurat_obj, reduction = name2D, group.by = "cluster_cell_type") +
Seurat::NoAxes() +
ggplot2::scale_color_manual(values = unlist(color_markers),
breaks = names(color_markers)) +
ggplot2::theme(aspect.ratio = 1,
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
clusters = Seurat::DimPlot(seurat_obj, group.by = "seurat_clusters", label = TRUE,
reduction = name2D) +
Seurat::NoAxes() + ggplot2::ggtitle("Clustering") +
ggplot2::theme(aspect.ratio = 1,
plot.title = element_text(hjust = 0.5))
clusters | cell_type
In this section, we’ll use the ISO_SWITCH_ALL function
from Isoswitch to analyze isoform switches between cell types. we refer
to “isoform switches.” when two isoforms of the same gene are markers of
different cell clusters. The ISO_SWITCH_ALL function uses
Seurat’s FindMarkers function to identify marker isoforms
for each cluster in a gene. It outputs a data frame, listing each
transcript marker in a separate row.
clusters = levels(seurat_obj@active.ident)
switch_markers = isoswitch::ISO_SWITCH_ALL(seurat_obj, clusters, assay = "ISO",
min.pct = 0, logfc.threshold = 0.40)
We save the list of DE isoforms :
saveRDS(switch_markers, file = paste0(out_dir, "./datasets/", "switch_markers_combined_all.rds"))
We load the DE isoforms table :
switch_markers = readRDS(paste0(out_dir, "/datasets/switch_markers_combined_all.rds"))
In this section, we explore the switch_marker
dataframe.
List of isoform switchs:
isoswitch::gene_switch_table(switch_markers)
Can you identify the genes and cell types that show the most significant changes in isoform expression?
To facilitate the graphical interpretation of the results:
plot_marker_matrix() produces a heatmap of number of
unique genes per contrast between clustersplot_marker_score() produces a volcano plot showing
p-values and average logFC for each gene with an isoform switchWe visualize the number of DE isoforms between each cluster :
pl1 = isoswitch::plot_marker_matrix(seurat_obj, switch_markers)
pl2 = isoswitch::plot_marker_score(seurat_obj, switch_markers)
pl1 | pl2
Individual volcano plots for each cluster.
isoswitch::plot_marker_score(seurat_obj, switch_markers, facet=TRUE, ncol=3)
In this section, we visualize isoforms expression level for specific genes.
First, we set the assay to “ISO” because all expression levels are stored there.
Seurat::DefaultAssay(seurat_obj) = "ISO"
Let’s take the example of Clathrin light chain A (Clta) gene isoform expression switch during neuronal maturation:
see_isoforms_expression(seurat_obj, "Clta", c(0,4)) %>%
patchwork::wrap_plots(., ncol = 3)
Seurat::FeaturePlot(seurat_obj, reduction = name2D, features = c("Clta..ENSMUST00000107849", "Clta..ENSMUST00000170241"), blend = TRUE, order = TRUE)
In this section, we generate a small report.
We load annotation files: The Isoswitch documentation provides a comprehensive guide on how to build your annotation metadata.
gtf_df <- readRDS(paste0(data_dir,"annotation/gtf_rds.rds"))
gene_metadata <- readRDS(paste0(data_dir,"annotation/gene_metadata.rds"))
transcript_metadata <- readRDS(paste0(data_dir,"annotation/transcript_metadata.rds"))
For Clta gene:
gene_switches <- isoswitch::compute_switches(switch_markers, gene="Clta")
knitr::kable(isoswitch::format_switch_table(gene_switches))
| geneId | t1 | c1 | p1 | log2fc1 | t2 | c2 | p2 | log2fc2 |
|---|---|---|---|---|---|---|---|---|
| Clta | ENSMUST00000107849 | Glutamatergic | 2.69e-23 | 5.35 | ENSMUST00000170241 | radial_glia | 4.11e-20 | -2.29 |
| Clta | ENSMUST00000107849 | Glutamatergic | 4.75e-27 | 2.35 | ENSMUST00000170241 | cyclin_radial_glia | 7.11e-20 | -2.03 |
| Clta | ENSMUST00000107849 | imature_Glutamatergic | 1.39e-20 | 5.03 | ENSMUST00000170241 | radial_glia | 3.87e-12 | -1.55 |
| Clta | ENSMUST00000107849 | GABAergic | 1.09e-18 | 5.54 | ENSMUST00000170241 | radial_glia | 1.56e-11 | -1.73 |
| Clta | ENSMUST00000107849 | GABAergic | 2.35e-19 | 2.54 | ENSMUST00000170241 | cyclin_radial_glia | 1.34e-10 | -1.46 |
| Clta | ENSMUST00000107849 | imature_Glutamatergic | 9.33e-20 | 2.03 | ENSMUST00000170241 | cyclin_radial_glia | 1.54e-10 | -1.28 |
| Clta | ENSMUST00000170241 | intermediate_progenitors | 3.60e-10 | -2.11 | ENSMUST00000107849 | Glutamatergic | 4.35e-09 | 1.21 |
| Clta | ENSMUST00000107849 | GABAergic | 3.37e-06 | 1.40 | ENSMUST00000170241 | intermediate_progenitors | 1.55e-05 | -1.54 |
| Clta | ENSMUST00000107849 | GABAergic | 1.09e-18 | 5.54 | ENSMUST00000107845 | radial_glia | 8.57e-05 | -2.16 |
| Clta | ENSMUST00000170241 | intermediate_progenitors | 1.50e-04 | -1.36 | ENSMUST00000107849 | imature_Glutamatergic | 8.12e-03 | 0.89 |
| Clta | ENSMUST00000107849 | imature_Glutamatergic | 1.39e-20 | 5.03 | ENSMUST00000107845 | radial_glia | 1.57e-04 | -1.84 |
| Clta | ENSMUST00000107849 | Glutamatergic | 2.69e-23 | 5.35 | ENSMUST00000107845 | radial_glia | 4.01e-03 | -1.94 |
| Clta | ENSMUST00000107849 | imature_GABAergic | 5.20e-07 | 4.49 | ENSMUST00000170241 | radial_glia | 7.68e-03 | -0.81 |
| Clta | ENSMUST00000107851 | Cajal_Retzius | 2.21e-02 | 6.09 | ENSMUST00000170241 | radial_glia | 2.58e-02 | -2.45 |
| Clta | ENSMUST00000107849 | Cajal_Retzius | 2.43e-04 | 4.38 | ENSMUST00000170241 | radial_glia | 2.58e-02 | -2.45 |
| Clta | ENSMUST00000107849 | Glutamatergic | 2.46e-09 | -0.85 | ENSMUST00000170241 | imature_GABAergic | 2.70e-02 | 1.49 |
| Clta | ENSMUST00000107849 | imature_GABAergic | 5.20e-07 | 4.49 | ENSMUST00000107845 | radial_glia | 2.97e-02 | -1.38 |
isoswitch::isoswitch_report(seurat_obj, "ISO", gene="Clta", marker_list=switch_markers, gtf_df, transcript_metadata)
## [1] "Locus"
## [1] "Junction"
## [1] "umi counts"
## [1] "dotplot"
## [1] "PATCHWORK"
We save the list of DE isoforms :
saveRDS(switch_markers, file = paste0(out_dir, "./datasets/", "switch_markers_combined_all.rds"))
## R version 4.2.2 (2022-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS 14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] reactable_0.4.4 ggplot2_3.4.0 patchwork_1.1.2 dplyr_1.1.4
##
## loaded via a namespace (and not attached):
## [1] backports_1.4.1 spam_2.10-0 Hmisc_4.7-2
## [4] plyr_1.8.8 igraph_2.0.1.1 lazyeval_0.2.2
## [7] sp_1.5-1 splines_4.2.2 RcppHNSW_0.5.0
## [10] crosstalk_1.2.0 listenv_0.9.0 scattermore_1.2
## [13] digest_0.6.31 htmltools_0.5.4 fansi_1.0.3
## [16] checkmate_2.1.0 magrittr_2.0.3 tensor_1.5
## [19] cluster_2.1.4 ROCR_1.0-11 limma_3.54.2
## [22] globals_0.16.2 matrixStats_0.63.0 spatstat.sparse_3.0-0
## [25] jpeg_0.1-10 colorspace_2.0-3 ggrepel_0.9.2
## [28] xfun_0.36 jsonlite_1.8.4 progressr_0.12.0
## [31] spatstat.data_3.0-0 survival_3.4-0 zoo_1.8-11
## [34] glue_1.6.2 polyclip_1.10-4 gtable_0.3.1
## [37] leiden_0.4.3 ggtranscript_0.99.9 future.apply_1.10.0
## [40] isoswitch_0.0.0.9000 abind_1.4-5 scales_1.2.1
## [43] DBI_1.1.3 spatstat.random_3.0-1 miniUI_0.1.1.1
## [46] Rcpp_1.0.9 viridisLite_0.4.2 xtable_1.8-4
## [49] htmlTable_2.4.1 reticulate_1.34.0 foreign_0.8-84
## [52] dotCall64_1.1-0 Formula_1.2-4 htmlwidgets_1.6.0
## [55] httr_1.4.4 RColorBrewer_1.1-3 ellipsis_0.3.2
## [58] Seurat_5.0.0 ica_1.0-3 pkgconfig_2.0.3
## [61] farver_2.1.1 nnet_7.3-18 sass_0.4.4
## [64] uwot_0.1.14 deldir_1.0-6 utf8_1.2.2
## [67] tidyselect_1.2.0 labeling_0.4.2 rlang_1.1.3
## [70] reshape2_1.4.4 later_1.3.0 reactR_0.5.0
## [73] munsell_0.5.0 tools_4.2.2 cachem_1.0.6
## [76] cli_3.6.2 generics_0.1.3 ggridges_0.5.4
## [79] evaluate_0.19 stringr_1.5.0 fastmap_1.1.0
## [82] yaml_2.3.6 goftest_1.2-3 knitr_1.41
## [85] fitdistrplus_1.1-8 purrr_1.0.2 RANN_2.6.1
## [88] pbapply_1.6-0 future_1.30.0 nlme_3.1-161
## [91] mime_0.12 compiler_4.2.2 rstudioapi_0.14
## [94] plotly_4.10.1 png_0.1-8 spatstat.utils_3.0-4
## [97] tibble_3.2.1 bslib_0.4.2 stringi_1.7.8
## [100] highr_0.10 RSpectra_0.16-1 forcats_0.5.2
## [103] lattice_0.20-45 Matrix_1.6-2 vctrs_0.6.5
## [106] pillar_1.9.0 lifecycle_1.0.3 spatstat.geom_3.0-3
## [109] lmtest_0.9-40 jquerylib_0.1.4 RcppAnnoy_0.0.20
## [112] data.table_1.14.6 cowplot_1.1.1 irlba_2.3.5.1
## [115] httpuv_1.6.7 R6_2.5.1 latticeExtra_0.6-30
## [118] promises_1.2.0.1 KernSmooth_2.23-20 gridExtra_2.3
## [121] parallelly_1.33.0 codetools_0.2-18 assertthat_0.2.1
## [124] fastDummies_1.7.3 MASS_7.3-58.1 withr_2.5.0
## [127] SeuratObject_5.0.0 sctransform_0.4.1 parallel_4.2.2
## [130] reactablefmtr_2.0.0 grid_4.2.2 rpart_4.1.19
## [133] tidyr_1.3.1 rmarkdown_2.19 Rtsne_0.16
## [136] spatstat.explore_3.0-5 shiny_1.7.4 base64enc_0.1-3
## [139] interp_1.1-3